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We discuss the analytical solution through the cavity method of a mean field model that displays 
at the same time an ideal glass transition and a set of jamming points. We establish the equations 
describing this system, and we discuss some approximate analytical solutions and a numerical strat- 
egy to solve them exactly. We compare these methods and we get insight into the reliability of the 
theory for the description of finite dimensional hard spheres. 
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FIG. 1: Schematic mean-field phase diagram of hard spheres in three dimensions, see the text and [5] for a detailed 
description. 
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I. INTRODUCTION 

The theoretical investigation of the glass transition and its relation to jamming in hard sphere systems 
has made considerable progress in the last 30 years [1-5]. This has been possible mainly because of the 
powerful analogy between jammed states and inherent structures [3, 6-8] and of the development of methods 
based on spin glass theory [9, 10] to describe the glass transition of particle systems. This progress led to 
the proposal that amorphous jammed states of hard spheres can be thought of as the states obtained in the 
infinite pressure limit of metastable glasses, and therefore described using tools of (metastable-)equilibrium 
statistical mechanics. 

The phase diagram of hard spheres that results from these mean- field studies is summarized in Fig. 1, 
where we plot the pressure as a function of the packing fraction ip which is the fraction of space covered 
by the spheres. The full black line represents the equilibrium phase diagram with the liquid-to-crystal 
transition. If this transition can be avoided (by compressing fast enough or by introducing some degree of 
polydispersity) , one enters into a metastable liquid phase. The nature of this metastable liquid changes at 
p = ifd- It consists of a single ergodic state for ip < ip^- When tp > <pa, the available phase space splits 
into many glassy states. If the system is stuck in one of these states and compressed, it follows one of the 
glass branches of the phase diagram, until its pressure eventually diverges at some packing fraction ipj which 
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depends on the state. At density ip^ a thermodynamic glass transition happens (in the sense of mean field 
spin glasses [11]) towards an ideal glass. The pressure of the latter diverges at (pocp- In the inset, the 
complexity, i.e. the logarithm of the number of glassy states, is plotted as function of the jamming density 
ipj\ this approach predicts that there exist jammed states in a finite interval of density ipj £ [<fth, 1 Pk\- The 
boxes show a schematic picture of the (3TV-dimensional, where TV is the number of particles) phase space of 
the system: black configurations are allowed by the hard-core constraint, white ones are forbidden. In the 
supercooled liquid phase the allowed configurations form a connected domain; however, on approaching ip^ 
the connections between different metastable regions become smaller and smaller. Above </?k, they disappear 
in the thermodynamic limit and glassy states are well defined. 

The above mean-field picture has been obtained by a succession of works which start from the studies 
of some categories of spin-glasses with so-called 'one step replica symmetry breaking', and have gradually 
matured into analytic approximation tools for the theory of hard spheres (see [5] and references therein). 
A very interesting model has been introduced recently by Mari, Krzakala and Kurchan [12]. It displays 
exactly the phase diagram presented in Fig. 1: it undergoes an equilibrium glass transition and it has an 
interval of densities where it shows all the phenomenology which is now associated to jamming, like marginal 
mechanical stability and the associated presence of anomalous soft modes in the vibrational spectrum [13- 
15]. The model has been studied numerically in [12] in order to show the existence of separate glass and 
jamming transitions and to clarify to some extent the relation between the two. 

This model is interesting in that it is in principle solvable: it can be investigated by mean of modern 
methods that have been developed in the context of mean field spin glasses, the replica method [16] and 
the cavity method [17]. This investigation is the purpose of the present paper, where we derive the cavity 
equations that describe the model and we present some approximated analytical solutions to them, along with 
a detailed numerical resolution. Since it will turn out that the exact solution requires quite heavy numerical 
calculations (heavier than a direct Monte Carlo study of the model, at least for a moderate number of 
particles, such as the one performed in [12]), one might wonder why this solution is interesting at all. There 
are at least two reasons why this study is interesting, in our opinion. The first is that Monte Carlo methods 
are not able to access the deep glassy phase or the densest part of the jammed phase: they are confined 
to explore the region close to ipd (at equilibrium) and (p t h (at jamming). Therefore if one wants to study, 
for instance, how the properties of the packings change when going from <p t h to ^gcp, the exact solution is 
needed. Moreover, we will show that the cavity method allows to derive simple analytical approximations 
to the true solution. Similar approximations have been used to study finite dimensional hard spheres [5]; 
their investigation in the controlled setting of the present 'solvable' model allows to assess their reliability. 
Finally, there are some generic structures in the correlations of jammed packings that one would like to 
explain analytically. Our work is a first step in this direction. 

This paper is meant to be read by specialists in the field, so we did not make much attempt to explain in 
details the basis of the method. Recent complete reviews of the physical problem [5, 18-20] as well as of the 
method we used [17, 21] exist, and the reader is assumed to be familiar with these concepts. 

II. DEFINITIONS 

The model that we study in this paper is a simple generalization of the one introduced in [12], defined as 
follows. We consider a "factor graph" , namely a bipartite graph made by two types of nodes: variables and 
boxes. Each variable is connected to z boxes and each box is connected to p variables. In a system with TV 
variables the number of boxes is Nz/p and the total number of links (i.e. variable-box connections) is TVz. 
We will consider an ensemble of 'random regular' factor graphs where each graph satisfying this requirement 
has the same probability. A crucial properties of this ensemble, that allows for the solution of the model, 
is that in the thermodynamic limit TV — > oo almost all graphs are locally tree-like, in a sense that can be 
defined precisely [17]. 

Each variable is a vector Xi £ [0, l] d with periodic boundary conditions, where d is the dimension and 
i = 1, • • • , TV. In the following we denote by \xi — xj\ = \J Y^=i(\ x i ~~ x> j Imod i) 2 the distance between 
Xi and its closest periodic image of Xj. If we call xi. x u x j) the characteristic function of the hard sphere 



FIG. 2: An illustration of the model for p — 6, z = 3 and N = 8. Each white square is a box, each black dot is a 
variable (sphere). Each box contains all the spheres connected to it by a link. The sphere inside one box must not 
overlap (note that for z = 1 one obtains N/p systems of p hard spheres) . 

constraint (with periodic boundary conditions), i.e. x( x i> x j) = 1 if \%i — Xj\ > D and otherwise, then each 
box a = 1, • • ■ , Nz/p imposes the condition 

hp 

i<j 

where xf are the variables connected to box a. The partition function of the model is 

f Nz/p 

Z= dxf-dxN J| X (a) • (2) 

J a=l 

A pictorial description of the model is the following (see Fig. 2). Each box can be thought of as a cubic 
region [0, l] d with periodic boundary conditions. Each variable node i = 1, . . . , N represents a "sphere" of 
diameter D and this sphere appears in position Xi in all the z boxes to which the node is connected. On the 
other hand, each box contains exactly p spheres. The constraint is that, for each box, the p spheres present 
in the box do not overlap. 

The model therefore differs from a standard hard sphere model, since each sphere interacts only with a finite 
subset of neighbors, and the topology of the interaction network is fixed by the random graph construction 
described above. This structure is such that the model becomes a mean field model and is therefore exactly 
solvable, at least in principle, as we will discuss in the following. It is worth to note, however, that there 
are two "formal" limits where one gets back the standard hard sphere model: in the case z = 1 the model 
reduces to N/p independent systems of p hard spheres each, while for p = 2 and z — N — 1 one gets back a 
single system of N hard spheres. Note also that in [12] only the version with p = 2 has been studied. 

Our investigations showed, however, that the model defined above undergoes a "crystallization" phe- 
nomenon at high density: the spheres tend to localize around a discrete set of positions inside the unit box. 
This has been avoided in [12] by introducing a small degree of polydispersity of the size of spheres. Here, in 
the analytical treatment of the model, we do not need to use this trick since we can impose directly that the 
solutions are translationally invariant, therefore discarding all crystalline phase of the model. In this way 
one effectively restricts to the amorphous phases, but one should keep in mind that these are metastable 
with respect to the crystal in the true model. Another possibility to remove the non-translationally invariant 
phase is to introduce local "random shifts": on each link we introduce a quenched variable s a i € [0, l] d , such 
that the corresponding particle appears in the corresponding box translated by s j. On a tree with open 
boundary conditions, this will not change the model since one can always perform a change of variable to 
remove the shifts. In presence of loops however, the random shifts will frustrate the periodic order. But 
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since the cavity solution is based on local recursions, the solutions describing the model with random shifts 
will be the same as the translationally invariant solutions of the model without random shifts. A similar 
situation occurs when studying an antiferromagnetic model on a random graph: local recursion relations 
allow both an antiferromagnetic and an amorphous ordering. The former is irrelevant on a random graph 
because long loops of odd length frustrate the antiferromagnetic order. The antiferromagnetic system thus 
behaves like the spin glass in which the sign of the couplings are quenched random variables. See Ref. [25] 
for a more detailed discussion in the context of a very similar model. 

We define Vd(R) the volume of a d-dimensional hypersphere of radius R; then V s = Vd(D/2) = 2~ d Vd(D) 
is the volume of one hard sphere (since the spheres have diameter D), and ip = pV s is the packing fraction, 
that represents the fraction of the unit box that is covered by the p interacting spheres. It is trivial to check 
that there are no configurations with ip > 1 . The parameter that controls the packing fraction is the diameter 
D since the box size is fixed; for this reason in the following we will use directly the sphere diameter D as 
control parameter and label the different transitions as D K , D GCP , D d , etc. 

For a system of p hard spheres in d dimensions, we define the following quantities: 



dxi ■ ■ ■dx p \Yx{x l ,x ) , 

i<j 

9p(x-y) = j-j (y2 5< y x ~ x i) 6 (y ~ x j)j = f dx3 ■ ■ - dx pX(x,y,x 3} - ■ ■ ,x p ) , 



such that Zp is the partition function of p hard spheres (apart from a pi), and g® is related to the usual pair 
correlation function [22] by 

9(r) = • (4) 

For the following discussion, it will be useful to define 



/n 



(5) 



which is the so called void space or cavity volume, namely the volume available to insert an additional sphere 
in a box given the positions of n other spheres, {x\, ■ ■ ■ , x n }. 



III. CAVITY EQUATIONS 

The cavity method has now become a standard method to solve statistical models defined on random 
graphs. We will not explain here the method and refer the reader to [17, 23]. Here we only write the 
equations for our specific case. 



A. Bethe free energy 

We define by di the set of boxes connected to variable i, and by da the set of variables connected to box 
a. On each link we define two fields: </ ? a->i( a; t) i s the probability density of the variable X{ when connected 
only to the box o; ipi-t a (%i) 1S the probability density of the same variable when connected to all the boxes 
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in its neighborhood but a. Both are normalized to 1 and they satisfy the equations: 



b£di\a 

(6) 



— I II dx ^^a{xj) \ x{a) 
a ^ iJ \jeda\i J 



which can derived from the stationarity of the Bethe entropy: 
S = 



l0g J dx i*Pi->a(x l )<fa^i(x l )+^2\0g J I Y[ dXjlpj^aiXj) J x(fl)+^log J dx l ]J fai) 
links a—i a \j£da J i aGdi 

' (7) 

These equations have the general form of the cavity (or Bethe) equations that can be derived for any model 
with local interactions [17]. With respect to previous studies of frustrated systems with the cavity method, 
the main difference here (and the main source of difficulty) is the fact that the variables x are continuous. 
Although the Bethe free energy is not variational in general, it has the property that the cavity equations 
can be obtained imposing its stationarity with respect to the cavity fields. In some special cases one can 
argue that it provides indeed an upper or lower bound to the true free energy, but a proof of this is still 
lacking. 



B. Replica symmetric cavity equations 

The replica symmetric (RS) equations for such a regular graph are trivially obtained by dropping the 
spatial dependence of the fields. In this case we use the notation = Z a ^i and Z^, = Zi^ ai and we get 



1 ft"- 1 \ (8) 

= ^Ws J \Y[ dx M x j) J X(x,xi,--- ,x p -i) , 
and the RS entropy per particle is 

/—..{../(IH - 

These equations admit the trivial translationally invariant solution tjj(x) — (p(x) — 1 with Z^ s = 1 and 

Zf = J dxj dx^ixj^j X (x,xi, ■ ■ ■ ,x P -i) = Zl , (10) 

that is the partition function of p Hard Spheres in the unit box. Therefore the entropy of the RS phase is 

S RS = -\ogZ° p . (11) 
P 
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C. 1-Step replica symmetry breaking cavity equations 



In the standard interpretation [17], the glass phase is signaled by the appearance of multiple solutions 
V'i^a' ^aXv °f ^q. (6). Each of these solutions represents a glass state with entropy s a given by the 
Bethc entropy (7) computed on the corresponding set of fields. Although one does not have direct access 
to individual glassy solutions (since the direct numerical solution of the Bethe equations by iteration on a 
single graph is extremely unstable in this region), a statistical treatment of the properties of the solutions 
in this regime exists and goes under the name of 1-step replica symmetry breaking (1RSB) description [23]. 
It is based on an entropy S(m) which is the sum over all solutions a of the corresponding partition function 
Z a = e NSa to power m [9] . The latter is computed by looking to the evolution of the solutions of the Bethe 
equations under an iteration that adds one more variable to the graph [23], or more simply by introducing 
an auxiliary model and assuming that a RS description holds for that model [17]. We do not discuss here 
these derivations and only report the resulting equations for our model, which are the following: 

S{m) = — logj^ Z™ = ms(m) + S(m) = -zS Unk (m) + -S box (m) + S site (m) 



Sunk(m) = log J dP[i>]dP[<p] J dxi/>(x)ip(x) 
S hox (m) = log J dP[^i] • • • dT[4> p ] J |n M*M 
S site (m) = log J dV[yi] ■ ■ ■ dV[ip z ] J dxf[<Pi( x ) 



P 

lGg«fc> 



The stationarity of this function with respect to V[ip] and V[(p] gives the 1RSB equations: 



i 

Zip 



z-1 



n win]* 
i=i 

p-i 

jjdvms 



ip(x) 



^-TT^(x) 



<p(x) 



Y~ J ndxjipj{x 3 ) X {x,xi, 



, Xp^lj 



where the normalization constants are 

Z^[tpi,--- ,<p z -i] 

Z v [ipi, ■ ■ ■ ,ip P -i] 
ZtI> = {{Zip)" 1 ) , 

z v = <(^r> • 



da; JJipi(a;) , 

i 

= J dxY\_dxjipj(xj)x{x,xi,- 



i Xp— i) , 



The internal entropy can then be written, using the standard method of [9], as 



s(m) 



dS(m) 
dm 



(Zunk^&Zlink) 



(z, 



link' 



zJZ, 
P 



bo JogZ box ) | (Z™ e logZ site ) 



box/ 



(12) 



(13) 



(14) 



(15) 



and the complexity is £(m) = S(m) — ms(m). The parameter m is the 1RSB parameter, whose equilibrium 
value must be fixed imposing that the replicated entropy is stationary [16]. 



IV. THE STABILITY OF THE RS SOLUTION 
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To study the stability of the RS phase we perturb around it: 

^ a (x) = l + Ae- ikx+i9 ^ , (16) 

and look at the linear stability of A assuming that the phase 9 is random, i.e. when substituting in the 
right hand side of (8) each ip get a random independent phase. This is done in order to enforce translational 
invariance, otherwise we would study the instability towards modulated phases, which is indeed interesting 
but we do not consider here, for reasons discussed in the introduction. Note that we have k — 2ir(m, • • • , rid), 
where m are integer numbers. Then at first order we have 



i< "' '" l^XX / \ ,x p )e- ikx >+»^" . (17) 

P o=l j=l J 



Now we can bring the factor e lkx on the other side and integrate over x; moreover we take the square and 
use that the 0j-> a are random and uncorrelated and we obtain the final result 



A 2 = A 2 (z-l)(p-l) 

Defining 



J dx!---dx P x{xi,--- ,x p -!)e 



2 

ik(xi — X2) 



(18) 



g° p (k) = j dxdye ik ^g p (x-y) = ^Jdx 1 ---dx pX (xu---,x p )e ik ^-^, (19) 



the stability condition is 



V(p-l)(z-l)\g p (k)\ < 1 , Vfc = 27r(ni,---,n < ,) 9 fe0. (20) 

Hence from the knowledge of Z° and g p (k) we can compute the RS entropy and the stability of the RS 
solution. 



A. Results for p = 2, any dimension 

For p = 2, k 7^ and D < 1/2, we have simply g%(x — y) = x(x — y)/(l — Vd{D)) and 



9°2(k) 



f e ik * X (x) _ f e ik *0(\x\<D) (2nD\ d/2 J d/2 (kD) 

A 0il] - l-V d {D) / 1/2 , 1/2]d l-V d {D) " \ k) l-V d (D)- { > 



One can show that for the values of D we are interested in, the maximum of g®{k) is assumed for k = 2n, 
i.e. the smallest k. Then the condition on D is 

D d/2 J d/ 2(2nD) _J_ 
l-V d {D) ~ Jz—\- 

In the limit z — > 00, as D is small, we can use J n {x) ~ (x/2) n /T(n + 1), and neglecting the denominator 

D^J^nD) ««/ 2 D« 1 

i-v^(D) ~fp+i) d( ) -W ( } 
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B. Results for d = 1, any p 

In d = 1 we get, from the exact solution 
Zl = [l-pDf- 1 , 

i p ~ 2 (24) 
9pW = — E e " W iF 1 [l + n-p--i{l- P D)k] . 

where iFi[a;b;z] is the confluent hypergeometric function of the first kind. Also in this case the lowest k 
becomes unstable in the first place. 



C. Results for d = 2 and p = 3 

As a last interesting case, we consider d = 2 and p — 3. In the following for simplicity we consider D < 1/4 
to avoid problems coming from periodic boundary conditions. 

We start by the computation of the partition function Z\ of three spheres in a box, which can be done 
using the standard virial expansion. For convenience we fix the first sphere, as well as the origin of the 
coordinate frame, in the center of the box. The center of the second sphere can be anywhere in the box 
outside a disk of radius D centered in the origin. Given the position of the second sphere, the third sphere 
can be anywhere outside the union of two disks centered around the first two spheres. 

If the second sphere is at distance r — |x 2 — x\\ from the origin x\ = 0, the free volume accessible to the 
third sphere is 



r r 



v 2 (x u x 2 ) = l-27TD 2 + 0(2D-r)D 2 ^2 arccos — - — ^ 4 - — j (25) 

This has to be integrated over the position of the second sphere. There are three possible cases: 

1. re [D, 2D]; in this case the first and second exclusion spheres have an overlap, and the second sphere 
can rotate at any angle without hitting the boundary of the box. Therefore one has 



p2D 

Z°(l) = 2ir dr 

J D 



1 - 2nD 2 + D 2 ( 2 arccos ^— - 



2D 2D 




(26) 



2. r£ [2D, 1/2] (recall that the box has side 1 so r is at most 1/2); in this case the first and second 
exclusion spheres have no overlap, and the second sphere can rotate at any angle, therefore 

,1/2 

Z°(2) = 2ir drr(l-2irD 2 ) (27) 

J2D 

3. r e [1/2,-^/2/2]; also in this case there is no overlap contribution, but the second sphere can only 
be at some angles because of the cubic shape of the box. The total angle that can be spanned is 
8(7r/4 — arccos(l/(2r))), therefore 

"V2/2 f / 1 

arccos — ) ) (28) 



Z§(3) = 8 / drr(l-2nD 2 

Jl/2 



4 V 2r 



All the integrals can be evaluated and summing the three contributions one gets the final result 

Z% = l-3irD 2 + jnD i (3V3 + 8n) , D < 1/4 . (29) 
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We also need the value of the pair correlation at contact, g®(D). Following the same reasoning this is 
given by 



0(D) = V2(r = D) = 1 *nU + L> y 3 2 J 



D < 1/4 . 



(30) 



Finally, 53(2; — y) = v 2 (x, y)/Z®, from which one can compute g®(k) numerically and determine the stability 
of the RS solution. 



THE GAUSSIAN APPROXIMATION 



We now introduce an approximation to describe the 1RSB phase of the model. We assume that the fields 
ipj (x) and ipi (x) are localized around a position which is randomly distributed in the box (this maintains the 
global translational invariance). This Ansatz, of course, is not a solution of the 1RSB equations. However, we 
expect that it provides a reasonable estimate of S(m), which is expected to become more and more accurate 
for large connectivity and close to the random close-packing point. Moreover, we will see in the following, 
that even if the variational nature of the replicated entropy cannot be proven, these approximations give 
upper bounds for Dk- For this reason we will refer from now on to these approximations as "variational" 
approximations. Note that if a variational approximation predicts that the Kauzmann radius is less than the 
radius where the RS solution is unstable, Dk < -Drs, then we know for sure that there is a discontinuous 
transition occuring at a value of D smaller than Drs. 

We assume a Gaussian shape for the fields, which leads to the following assumption for their distribution: 



PM = J 

m = J 



dXS 



dX5 



(2irA) d / 2 



ip(x) - 



(31) 



{2n5A) d / 2 

We substitute this Ansatz in the Bethc free energy (12) and determine the variational parameters A and 5 by 

its extremization. In the following we will use the definition 7a(x) = ^Ay 1 / 2 • Substituting the expressions 
above in (12), we obtain the following results: 



Slink 
Ssite 



: lOg 

log 



m- d / 2 [27r(l + 6)A] d < 1 - m V 2 



Jl-z)d/2 (l- m )d/2 



(2tt5A) 



-(l-m)(l-z)d/2 



(32) 



Note that Sbox does not depend on 8. Therefore we first write the contribution of Su n k and S S i te and optimize 
with respect to 5: 



S site - zSunk = ~(1 - m)log(27ivl) + ^ log to + ^(1 - to) log 



(i + sy- 



(33) 



The optimization is straightforward and gives 5 — z — 1 as expected from the first Eq. (6). The optimized 
result is 



|(1 - to) log(27r^) + ^ log to + |(1 - m)(z - 1) log 



The last term to be computed is Sbox, which has the form: 



i-I 

z 



Sbox = log / dXt ■ ■ ■ dX v 



dxi ■ ■ •dx P 7 j4 (xi - Xi) ■ ■■'ya{xp - X p )x(xi, ■■ ■ ,x p ) 
Unfortunately this cannot be computed exactly and we have to resort to further approximations. 



(34) 



(35) 
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A. Small cage expansion, first order 

The small cage expansion proceeds as follows [5]. First we assume that m is an integer and write Sbox as: 

hp 



J box 



= log / dxx ■ ■ ■ dx p p(x l ) ■ ■ ■ p(x p )Y[x( x i,Xj) , 



(36) 



where x = (x\, ■ ■ ■ ,x m ) is the coordinate of a "molecule" made of m particles, xi x iV) = Il^Li xi x a,y a ), 
and p(x) = J dX ]\™ =1 1a ( x a — X). Observing that J dx 2 ■ ■ ■ dx m p(x) = 1, we write 



S box = log / dxi--- dx p p{x 1 ) ■ ■ ■ p(x p ) Y[[x{xi,Xj) - x{xu,xij) + x( x u,xij)] 

J i<j 



log 



hP 



hp 



hP 



J dx n ■ ■ ■dx lp Y[x(xii,xij) + J dx u ' ' ' dx lp | J| x( x w, xij>) J Q{xu - xy) 



i<j i<j 

where we omitted the second order in the development in series of x ~ Xi an d we defined 



(37) 



Q( x - V) = J dx x ■ ■ ■ dx m dyi ■ ■ ■ dy m p(x)p(y) 



| I x(x a ,y a ) - i 



a=2 



(38) 



In [5] it is shown that the second order gives a contribution O(A) and that at lowest order (see Appendix 
C3 of [5]) Q(r) = 2^/AQ a (m)5(r — D), where Qo{m) is a function of m defined in [5] as: 



Qo{m) = r [Q(t) m - 0(t)] ; 0(t) = hi + erf(i)] = -h f dxe 



(39) 



We get then 



^-logZO + fcH 



log^° 



p(p - 1) 2dVA 
2 IT" 



dxdyQ(x - y)gl{x - y) 



Q (m)g° p (D)V d (D) 



(40) 



and collecting all the terms we get 



S(m) =^(m - 1) log(2^A) + | logm + |(1 - m)(z - 1) log 



1-1 

z 



z z(p- 1) 2dVl 



(41) 



+ -logZ£ + 



Optimization with respect to A gives 



Qo(m)g°p(D)V d (D) . 



'A* = D 



1 — m 



Qo(m)z(p-l)V d (D)gO(D) ' 



and 



S(m) =^{m- l)log(27rA*) + ^logm + d(l - m) + ^(1 -m)(z- l)log 



(42) 



+ -logZ°. (43) 
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FIG. 3: Special values of the sphere radius as functions of z at p — 2 for different values of d in the Gaussian 
approximation: Drs beyond which the RS solution becomes unstable, Dgcp where the pressure diverges, and Dk 
where the Kauzmann transition takes place. When Dk < Drs the transition is necessarily first order. 



In particular, using the results Qo(m — > 0) ~ y/n/4m and Qo(m ~ 1) = Qo X (1 — m) with Qo = 0.638 [5], 
one can show that this expression trivially reduces to the RS entropy (11) for m = 1, and that 



Sj = lim S(m) — —dlog 

m— >0 



2V2D 



S e o = — lim m d m [S(m)/m] 

TCI— ^1 



z(p-l)V d (D)g°(D) 

dlog 



+ d+-(z-l)log 



1-i 

z 



d, 2tt 
— log — 

2 6 e 



D 



z(p-l)V d (D)g°(D)Q 



-\ogZl 



+ -(z_l)log 



B. Results for p = 2, any dimension 



For p = 2 we have trivially Z$ = I - V d (D) and g°(x,y) = x( x ,v)/ Z 2> therefore g%{D) = 1/Z$. We get 



S(m) = ^(m — 1) log 



27rL> 2 (l-T/ d (L>)) 2 (1-m) 2 



+ d(l - m) + -(1 - m)(z - 1) log 



Mm) 2 
1" 



• log m 



■ io g [i-y rf (£»)] 



(44) 
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FIG. 4: Drs, Dgcp and Dk as functions of z for different values of p at d = 1 in the Gaussian approximation. 



and 



Ej = lim S(m) = ——log 



m— >-0 



8D 2 (1 - V d (D)f 



z 2 V d {Df 



+ d+-(z-l)log 



\og[l-V d {D)\ , 



lim m 2 d Tn [S (m) / m] — log 



m— >1 



2ttD 2 (1 - V d {D)f 
z^V d {DYQl 



d d, 
+ 2 + 2^- 1 ) 1 °S 



1 

1- - 

z 



+ -log[l-V d (D)} 



and Dk is defined by Y, eq = while -Dgcp is defined by Xj = 0. The results are reported in Fig. 3. 

C. Results for d = 1, any p 

Also in d = 1 the integrations can be performed for all p. We get 



9° P (D) = 



1 

1-pD 



(45) 
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Then 



and 



n(l-pD) 2 (1-m) 2 



4z 2 (p-l) 2 Qo(m) 2 



+ 2 lo S TO 



+ (1 - m) + -(1 - m)(z - 1) log 



;l0g 



■log 



2(1 -pDf 



z 2 {p-lf 
n(l-pD) 



1 

1 - - 

z 



1 - 



+ Z{P ^ log(l-pJ) 



z(p - 1) 



Iog(l -pD) , 



2 1 



2z 2 (p-l) 2 Q 2 



5 + ^-1) log 



P 



(46) 



log(l -pD) . 



The results are reported in Fig. 4. 



VI. THE DELTA APPROXIMATION 



In this section we introduce another variational approximation scheme, that we shall call the "delta 
approximation" . The motivation is that within the Gaussian Ansatz, A — > at jamming: therefore, both 
ip{x) and ip(x) become delta functions in this limit. We would therefore like to compute the free energy 
directly for delta function fields; we expect this to give a simpler expression of the free energy, that should 
be good close to jamming. The problem is that the Gaussian expressions are divergent for A — > unless m 
also goes to zero proportionally to A. This is due to the fact that both fields ip(x) and <p(x) become delta 
functions for A — > 0. We therefore construct here a different approximation by eliminating the field ip(x) 
and making a delta function Ansatz only for the field ip(x): in this way the field tp(x) is computed exactly 
and in particular it is not a delta function. 

One can show in general that by using equations (13), one can eliminate the field (f(x) and the replicated 
entropy can be equivalently written as 

S(m) - S slte , - Z -^—±±S box (47) 
where Sbox is defined as in Eq. (12) and 

S si te> =log J dV[ti]---dV[r p -i] J dx f[J dx k 1 ---dx k p _ 1 ^(x k 1 )---^_ 1 (x k p _ 1 ) X (x,x k 1 ,--- ,4^) 

(48) 

The "delta approximation" is then based on the following Ansatz for V(ip): 

V[^} = J dXS [iP(x) - S(x - X)] , (49) 

namely on each site i the probability of the variable Xi is a delta function centered in a i.i.d. random point. 
Under approximation (49), the replicated entropy becomes 

(z \ m 

J dx n X(x, X k , ■ ■ ■ , X k _M - z{P p 1} log J dX! ■ ■ ■ dX pX (X u ■■■ ,x p ) 

= lo s / (il dx " ■ ■ ■ dx p-ix(xt ■ ■ ■ , x k p _^ f z (p_i) (x\ ■ ■ ■ x;_,r z -^-^- log z° p , 

(50) 
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recalling the definition of v n in Eq. (5). Introducing the normalized measure of n spheres in a unit box, 

, , dxi ■■■dx n x(xi ■■■x n ) 

d/j,(xi ■ ■ ■ x n ) = ^ , (51) 

we can rewrite S(m) given in Eq. (50) in the equivalent form 

S(m) = log J (f[ d^X* ■ ■ ■ X k p _^ [v^-v (Xl---, X;_j\ m + z log Z° p _, - Z -^-^ log Z° p . (52) 

In the following we study this expression for several specific values of p and d. In this section we will derive 
the expressions for the complexity, and in section VIII we will present the results together with a comparison 
with numerical resolution of the cavity equations. Note that for m — 1 one can easily show that S(m) 
given above is equal to the RS entropy (11), which is an important requirement for the consistency of this 
approximation. 



A. One dimension 

1. Results for p — 2 

We first consider the simplest case, namely one spatial dimension and only two-particles-in-a-box interac- 
tions (p = 2). Since Z\ = 1 and Z% = (1 - 2D), we get 

S(m) = log f f[ dX t [v^X, ■ ■ ■ X z )] m - Z - log(l - 2D) . (53) 
J i=i " 

We have therefore to compute the probability distribution P z (v) of the void space left in [0, 1] for the insertion 
of a new particle, after having put z particles in random positions {Xi}. Then we have 

S(m) = log / dv P z (v) v m - - log(l - 2D) . (54) 
Jo 2 

Note that v ranges from (no void space) to 1 — 2D (in the limiting case where all points Xi coincide) , and 
we expect that P z (v) = poS(v) + P z e9 (v) since a finite fraction of configurations have zero void space at large 
enough D. Since the delta function does not contribute to S(m), we will omit it from now on. 

In order to estimate P z {v) we can make the assumption that whenever v > 0, there is only one hole large 
enough to contribute to v (i.e. a hole whose length is bigger than 2D). The function P z (v) can then be 
easily evaluated in the following way. The hole that contributes to v must have length 2D + v, and must be 
delimited by two particles that we can choose in z(z — 1) different ways, since particles are distinguishable. 
We can put the first particle in x\ = and the second in x 2 = 2D + v (integration over x\ can be omitted 
since it gives a factor of 1 , the length of the box) . The remaining z — 2 particles must be in the space between 
X2 and 1, therefore giving a contribution (1 — 2D — v) z ~ 2 . Therefore, within the one-hole approximation, 
we get P z (v) — z(z — 1)(1 — 2D — v) z ~ 2 . We notice that the total probability of v > must be smaller then 
one since some configurations might have v = 0. This gives the condition 

1-2Z5 

dvP z (v) =z(l-2D) z - 1 < 1 => D > (1 -z- 1 /^- 1 ))/2 , (55) 

which gives an estimate of the limits of validity of the one-hole approximation. 

Plugging the result for P z (v) in Eq. (54), we get an approximate formula for the replicated free energy 
which depends on z and D, 

5( '»> = '° g ( r( ° ?('+', 7" ) + ('» - 1 + 1) log<1 - 2D » • (56 » 



L 



Recall that E 
get 



eg 



[m 2 d m (S(m)/m) 
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and that Dk is the point where the latter quantity vanishes. We 



Se, =2^ + ^^(1-21)) 

q=2 q 



D K = 



1 — e "~ 2 ^ 



(57) 



On the other hand, T,j = S(m = 0) and it vanishes at the close packing diameter Dgcp- We get 



- log(z) + i-i log(l - 2D) , D GC p = 1 



1 _ z -V("-2) 



(58) 



The complexity curve can be obtained explicitely, using £ = — m 2 d rn (S(m)/m) and s = d m S(m), which 
gives the parametric representation: 



a = log(l - 2D) - ^ , 
^— f m + q 

9=1 



z - 2 

E = — — log(l - 2D) + log 



r(z + i)r(m + i) 

T(z + m) 



(59) 



+ m E 



^— ' m + g 

9=1 



One can check easily that both critical diameters Dk and Dgcp are well within the region of validity of the 
one-hole approximation given by Eq. (55), and they scale as Dk,Dgcp ~ logz/z in the large connectivity 
limit. The values of Dk and -Dgcp can De compared to the stability of the RS solution (which scales as 

D s ~ 1/y/z). 



2. Results for p — 3 

We now consider the three-particles- in-a-box case p = 3, still for d = 1. Since Z\ = 1 — 2D and 
Z§ = (1 - 3D) 2 , we get from Eq. (52): 

S(m) = log / dv P 2 , z (v) v m +z log(l - 2D) - — log(l - 3D) . (60) 
Jo 

where now P2, z {v) is the probability distribution of the void space in [0, 1] for the insertion of a new particle, 
after having thrown at random z pairs of particles, each pair being at distance bigger than D. The latter 
ranges from (no void space) to 1 — 3D (in the case where each pair is exactly at distance D and superposed 
to all the others). 

Within the same one-hole approximation, we can approximate P2,z{v) as follows. The hole must have 
length L = 2D + v. We have to distinguish between two different situations: i) The hole is made by the same 
couple of particle; ii) The hole is made by two different couples. In the case i) we have z ways of choosing 
the couple. We fix then one of the two particles of the couple in and the other one in L (which gives an 
extra factor 2). Finally the other z — 1 couples of particles must be in the interval [L, 1] with the conditions 

that they are pairwise compatible, which gives a factor f(L,D) = J^dx dy \( x ,y) = (1 — L — D) 2 for 
each pair. With this definition the contribution due to the same couple finally reads: 2zf(L,D) z ~ 1 . In the 
case ii), instead, we can fix one particle of one couples in (we have 2z ways to choose it) and one particle 
of another couple in L (we have 2{z — 1) ways of choosing it). The free particle of the first couple must be 
in [L, 1 — D], due to the condition that it is compatible with its partner which has been fixed in 0. This 
gives a contribution (1 — L — D). An analogous contribution comes from the the free particle of the second 
couple, which must be in the interval [L + D, 1]. The other z — 2 couples must be in the interval [L, 1] and 
must satisfy the compatibility condition, and therefore give a contribution f(L,D) z ~ 2 . The sum of the two 



17 

contributions is (4z 2 — 2z)(l — L — D) 2 ( z ~ x \ and it has to be normalized by the total integral (1 — 2D) Z ; 
going back to v = L — 2D we get 



^2,» 



2z{2z- 1)(1 -v -3D) 



2(z-l) 



(1 - 2D) 2 



As in the previous case we get the condition 

r-l-3D 



Jo 



(1 - 2D) Z 

which gives a lower limit of validity in D of the one-hole approximation. 
Plugging this results in Eq. (60) we get for the replicated entropy 



(61) 



(62) 



S(m) = log 



r(™ + i)r(i + 2z) 



T(m + 2z) 



m-l-y ) logo -- m 



from which we get 



and 



q=2 q 



£, = log(2z) + HlJ! l g(l - 3D) 



l r 



D K = 



1 — g 22-3 2^ij = 2 q 



1 - (2z) 



-3/(22-3) 



(63) 



(64) 



(65) 



We checked that both Dqcp and Dk are well within the region of validity of the one-hole approximation; 
actually, the value of the left hand side of Eq. (62) never exceeds 0.1. Again, Dqcp an d Dk are found to 
scale as 2 log z/ z for large z. 

3. Conjecture for arbitrary p (2, 3, • • • ,oo) 
A comparison of Eqs. (57) and (63) and of Eqs. (58) and (64) allows to guess the form for general p: 



(p-!) z 1 / , \ 
\ q = E - + ip - 1)Z - P lo g (l- P B), 



Iog((p - l)z) + 



(p- l)z~: 



log(l - pD) 



f K = — 

P 



1 — e (p-i)z-p 
1 

p — - 



v^(p-l)« 1 
2^9=2 , 



(66) 



1 - ((p- 1) Z )-P/((P-!) Z -P) 



however we did not attempt to provide a proof of this conjecture. 



B. Two dimensions 



In the d = 2 case we cannot compute S(m) analytically and we must resort to a numerical evaluation. 
The numerical algorithm consists in writing a routine that is able to compute the void space v n , defined in 
Eq. (5), left by n disks centered in a set of positions {X}. We used an adaptation of the algorithm described 
in [24] that works as follows: 



We start by a grid of squares of side A -c D (typically A = 1/100). These squares are considered as 
particular cases of convex polygons. 
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• We add disks X\ ■ ■ ■ X n sequentially. 

• Each time a disk is added, we check if a given polygon is entirely contained in the disk. In this case it 
is removed from the grid. 

• Next we consider the polygons that intersect the boundary of the new disk. We approximate the 
boundary of the void space left in the old polygon by a new polygon, by approximating the boundary 
of the disk by a straight line (which is reasonable if A <C D, with error 0(A/ D) 2 ). The new polygon 
replaces the old one in the grid. 

• This construction is iterated until all disks have been placed. The area of the polygons that survived 
is computed easily using Eq. (1) of Ref. [24], and it gives the void space v n . 

The void space has to be averaged over the distribution Y\i=i d[i(X{ ■ ■ ■ dXp^), hence we must sample a 
configuration of p — 1 spheres in a box (and do this z times indepentently). This can be easily done for p = 2 
(one sphere, flat distribution) and p — 3 (put one sphere in the centre of the box, draw a second sphere 
outside it, then translate randomly both spheres). 

A correct sampling gives access to the void space distribution -P(w), that has the form P(v) = p 5(v) + 
P reg (v), as in one dimension. In the following we omit the delta term and only consider P reg (v), which 
therefore is not normalized to one (its integral gives the probability that v > 0). From this we can compute 
Eq. (52) as we did in one dimension: 

S(m) = log J dv P{v) v m + z log - Z{jP ~ ^ log Z° . (67) 

Similarly we get, using the relation J dvP(v)v = (v) = (Zp/Zp_ 1 ) z (which can be easily checked and also 
serves as a check of the correct sampling of P(v)), 

E e9 - - log Z° p - ( / dv P(v) vlogv, 

P \ Z p J J (68) 

- log J dv P(v) + z log Z p °_! - z{P p 1} log Z° p . 

Therefore both E eg and £j can be computed directly from P(v); from them we can determine the transition 
points £>k and Dqcp- 



VII. NUMERICAL SOLUTION OF THE EQUATIONS 



In the previous sections we described two analytical approximate methods yielding the phase diagram of 
the model. Beyond these analytical approaches, one can also develop some algorithms to solve the functional 
self-consistent 1RSB equations numerically. In this section we explain how it is possible to implement a 
numerical procedure to solve Eqs. (13) in the 1RSB phase for each value of the connectivities, z and p, of the 
diameter D, of the 1RSB parameter m and, in principle, of the spatial dimension d (in practice, numerical 
solutions can only be achieved in one and two dimensions). In order to do that we need representations of 
the cavity fields <p(x) and ip(x), and of the distributions V[<p] and V[tp], which can be treated by a computer. 

As far as the cavity fields are concerned, the simplest possibility is to discretize the volume [0, l] d where 
the functions (p(x) and ^(x) arc defined using a regular hyper-cubic grid with q bins per side of size 1/q. 
For instance, in one dimension we discretize the interval [0, 1] in q slices of length 1/q, and in two dimension 
we discretize the square box on a square lattice of q x q points. 

The coordinate in the box can assume a discrete set of values, i/q, with i being a rf-dimensional vector 
whose components arc integers between and q — 1, identifying the coordinate of the position of the center 
of the sphere in the box. If the position of the center of the sphere occupies a given site of the grid i, then 
all other sites of the lattice that are at Euclidean distance from i smaller than the diameter of the sphere 
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diameter as D c g = | 



D cannot be occupied by the center of another sphere (we call this number no)- The volume of the sphere 
in the discretized version of the model can be estimated as V s — ri£>/(2g) d , and the packing fraction as 
(p = pV s = pri£>/(2q) d . Since in the continuum limit V s — Vd(l)(D/2) d , we can then define an effective 

l/d 

. Note that in general Z? c ff 7^ D, and we take D e g as representative of the 
sphere diameter in the continuum limit. In particular, by symmetry, in d = 1 the number of excluded sites 
always has the form nn = 1 + 2a for integer a, and one has 

(«) 

In d = 2 the parameter no depends in an irregular manner on the choice of D (since the square lattice we 
use breaks the spherical symmetry) and one has in general 

(70) 

In the discretized version, the fields <p{x) and ^{x) are vectors of q d components (such that the sum of all 
components is equal to one), and the cavity equations, Eqs. (6), become a set of coupled algebraic equations 
for the q d components of the cavity fields, which can be easily solved numerically (of course, the numerical 
complexity of this step grows linearily with the number of components of the cavity fields, q d ). 

Note that the discretized version of the model is a generalization of a very important optimization problem 
known as the "random graph coloring" problem, where the number of colors corresponds to the number of 
components of the cavity fields q d . In particular, for no = and p = 2we recover the standard q-coloring 
problem, which has been deeply studied in the past few years, and whose properties and phase diagram are 
known in great details [25]. 

The continuum limit of the model is, of course, recovered for q — > 00. As a consequence, in order to make 
sure that the numerical results are reliable and that they are not affected by the discretization, we solve 
numerically the 1RSB equations using several values of q, and analyze the scaling properties of the numerical 
solutions with the number of bins. Moreover, one should note that for d > 1, partitioning the box using an 
hyper-cubic grid breaks the spherical symmetry down to some discrete symmetry. This makes the scaling 
towards the continuum limit in two dimensions more problematic than in one dimension (also because, due 
to the fact that the complexity of the numerical algorithm grows as q d , we are limited to smaller values of q 
ford = 2). 

Other numerical representations of the cavity fields were also possible. For instance, as f(x) and tp(x) are 
periodic functions in the interval [0, l] d , we could have performed a Fourier transformation of the recurrence 
equations keeping all the components up to a certain momentum, yielding a finite set of coupled algebraic 
equations for the Fourier coefficients of the cavity fields (similarily to what we did in Sec. IV to study the RS 
stability). However, it turns out that this strategy is not efficient in the most interesting region of the phase 
diagram, namely at high packing fraction where a 1RSB glass transition is found. Indeed here the cavity 
fields becomes extremely peaked (this is also the reason why the Gaussian and the delta approximation work 
very well), and the momentum cut-off needed to get accurate results becomes too big to be handled. 

Another possibility we could have employed, is to represent the fields as a population of delta functions, 
e.g. <p(x) = ^2 a c a 5(x — x a ). This strategy, which has the advantage that one does not need to discretize 
the space, has, on the other hand, the disadvantage that at each step of the iterative procedure, in order 
to generate a new field, one has to sample uniformly one point in the free space available for the insertion 
of a new particle, given the position of z(p — 1) neighboring particles in the box. This is trivial in d = 1, 
however in that case the discretized procedure work already well enough. In d = 2, this could be done using 
the algorithm described in Sec. VI B. However this algorithm is too slow to be used efficiently to this scope. 
Therefore in the following we will not explore further this representation. 
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A. The population dynamics algorithm 



Now, once that we dispose of the discretized representation of the cavity fields, we need to be able to 
implement a computational strategy to solve the 1RSB functional self-consistent equations, Eqs. (13), for 
any value of the connectivities, z and p, of the diameter of the spheres, D, and of the 1RSB parameter m. 
This step is quite standard in the context of the cavity method, and goes under the name of "population 
dynamics algorithm" [23]. The idea is to represent the probability distributions V[ip] and P[ip] as populations 
of M. representative cavity fields with some weights: 

M M 

P[<p] = Y,z25[<p{x)-<p a {x)], and Pty] = £ z% S^x) - i> a {x)) (71) 

a—l a—1 

As previously discussed, we need to consider only translationally invariant solution of Eqs. (13) in order to 
describe the glassy phase. A solution P[ip(x)] is translationally invariant if the property P[i[)(x+s)] = V[i[)(x)] 
holds for any s E [0, l] d , where ip(x + s) is an arbitrary translation (taking into account periodic boundary 
conditions) of ip(x). Since we represent the probability distribution P[ip] by a set of representative samples 
ip a (x), it is very easy to implement translational invariance. In principle, we would like to impose that 
if ip a (x) is one of the samples, then any translation of it is also contained in the set of samples with the 
same weight. But this is just equivalent to do the following: at each time we use a given sample i[)(x) as 
a representative of P[tp], we apply to it a "random shift", namely we extract a vector s uniformly in [0, l] d 
and we translate i/j(x) by s. In this way we impose translational invariance by hand. 
The population dynamics algorithm works in the following way: 

1) Pick at random p—1 fields ipi from the population P [ip] , according to their weights z^. Apply a random 
shift with fiat probability in [0, l] d to each of the cavity fields. 

2) Using Eq. (6), compute the new cavity field <p, along with its weight z v , which is given by the nor- 
malization in Eq. (14) to the power m, according to Eq. (13). Note that at high density, in the 1RSB 
phase, the cavity fields becomes extremely peaked. This implies that there exist some configurations of 
the p—1 fields ipi for which the new field ip is zero everywhere in [0, l] d . In this case the corresponding 
weight is zero and we have to reject it and restart the procedure. These events, which can cause a 
major slowing down of the algorithm, are called "rejection events" . 

3) Repeat 1) and 2) M. times, until a whole new population P new [<*p] is generated, and replace the old 
population with the new one (this kind of update is called in the context of population dynamics 
algorithm "parallel update"). 

4) Apply steps 1), 2), and 3) using the population P[tp] to generate a new V n ew[i>]- 

5) Repeat steps 1), 2), 3), and 4) until convergence, namely until the populations P[ip] and P[<p] are 
stationary. 

Once this process has converged, we can compute the average values of the link, the site and the box 
contribution to the 1RSB entropy, Eq. (12), from which one can obtain the complexity £(m). This allows 
to determine the equilibrium value of m* inside the 1RSB glassy phase as the point where S(m) has a 
minimum [9]. In practice, instead of computing the replicated entropy using Eq. (12), we can use another 
and equivalent formula (derived below) which is more advantageous from a numerical point of view. Indeed, 
using Eqs. (6) we can easily obtain the following relations (we omit the arguments of the functions Z): 

Zlink = = ■ (72) 
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Using these and Eqs. (13), one can rewrite the total and internal entropy as 

S(m) = (l-z + ^j S Hnk + Z -S v + SV 

s(m) = ( 1 - z + - ) si ink + -s v + (73) 



The computation of S v = \og(Z™) and 5^ = \og{Z™) is numerically less involved than S S i te and Sbox 
appearing in Eq. (12). Moreover, these contributions can be evaluated ondine during steps l)-5) of the 
population dynamics algorithm described above (we have just to compute the average value of Z™ and Z™ 
over all the M. attempts of generating a new cavity field), without requiring the implementation of any 
further step. 

Of course, representing the distributions V [ip] and V [ip] as populations of M. elements is an approximation 
which becomes exact only in the Ai — > oo limit. On the other hand, the numerical complexity of the 
population dynamics algorithm grows lincarily with M. In practice on has to find a good compromise 
between a value of M. small enough such that the execution time of the code stays reasonable, but big 
enough to avoid systematic corrections due to the finite size of the populations. In the present case, we find 
that M — 2 16 is close to the optimal value. 

Although we have produced a working version of the algorithm described above at any finite value of the 
1RSB parameter to, it turned out that the execution time is too big to get accurate results in a reasonable 
time. However, there are two special limits, namely m — > 1 and m — > 0, which describe respectively the 
physics at the Kauzmann point and in the close packing regime, where some semplifications arise which 
allow to perform the numerical study of the model in a more efficient way. These two limits are discussed 
below. 



B. Reconstruction: the limit m = 1 



In this section we consider the numerical solution of the 1RSB equations for to = 1. Recall that S(m = 1) 
gives back the equilibrium RS entropy of the system between the dynamical transition (where a non-RS 
solution of the 1RSB equations appears for the first time due to the emergence of glassy metastable states) 
and the Kauzmann point. In this limit, using the approach introduced in [26] which goes under the name 
of reconstruction method, also applied in a similar context to the coloring optimization problem in [25], the 
self-consistenf 1RSB equations can be simplified. Similarily to [25, 26], one can indeed introduce two new 
families of distributions over the cavity fields for each value of the variable x, defined as 



K x [ip] = i/>(x)V[il>] and TZ X 



(74) 



Using the previous definitions, the 1RSB cavity equations, Eqs. (13) can be rewritten in terms of these new 
distributions. Furthermore, imposing the translational invariance which implies that lZ x [ip(y)] = lZo[<p(y— x)] 
for all x we obtain the the self-consistent recursion relation for the new distributions which read: 



.2-1 r 1 

ftoM = ndfco[<pi]6 il>(x) - — Y[<pi(x) 

J i=i L ^ i 

Tl [ip} = j dfj,(xi---x p -i\0)Y[dRo[ipi}S <p{y)-— j J[ dyji}j{yj - Xj)x(y,yi, ■ ■ ■ , V P -\) 



(75) 



where 



dfi(xi 



Cp-l 



10) 



X(0,X!,--- ,Xp-i)dx 1 ---dx p -i 



(76) 
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From a numerical point of view, these latter equations are much easier to solve than Eqs. (13) for two reasons. 
First, no reweighting factor is present, which prevent the population to concentrate on few cavity fields with 
large weight. Second, rejection events cannot occur in this case. Indeed, for example, the procedure to 
generate a new field ip amounts to: 

1) Pick at random p— 1 fields ipi from the population 1Zo[ip]. Note that all the fields have the same weight 
in this representation. 

2) Pick p — 1 variables Xi,--- , x p _i in the interval [0, l] d satisfing the hard-sphere constraint 
X(0, xi, • • • , x p _i) with a fiat measure. 

3) Shift each of the p—1 chosen cavity fields by Xi. 

4) Using Eq. (6), compute the new cavity fields ip (again, note that there is no reweighting in this case), 
and insert the new field randomly into the population 7\LoM (this kind of update is called "serial 
update" and ensures a better convergence than the parallel one). 

Once the populations T^oM and T^-obP] have attained stationarity, we can compute the complexity of the 
system. Since the replicated entropy S(m = 1) equals the RS one, the complexity at m = 1 is given by 
S e «j = Sb,s — s ( m = !)• The internal entropy can be evaluated using Eqs. (15) and (13), where 

(Z Hnk log Zi in k) = J dn [^}dn [(p} log J dy ip(y)ip(y) 

(Z lp \ogZ lp ) = J JJdftofo] log J dyY[p(y) (77) 

i—1 "' i 

(Z^logZ^) = d^(xi---x p ^i\0)ZpY\dTZ a [^i}\og dyY[dyii/)i(yi-Xi)x{y,yi,---,y P -i)- 
•> »=i •* i 

From the complexity we can determine the Kauzmann point, which corresponds to the value Dk where S eg 
vanishes. 

In principle this method would also allow to determine the location of the dynamical transition, which is 
the first point where a non-RS solution of the 1RSB equations appear at m = 1. 

The results at m — 1 obtained with the reconstruction method will be discussed in Sec. VIII, and compared 
with the analytical approximations. 



C. Hard fields: the limit m — 



Also this specific limit yields a simplification of the numerical algorithm. The m — > limit corresponds in 
this context to the "close packing limit", since an inspection of the expression of the internal entropy s(m) 
shows that it goes to — oo as log(m), and the pressure diverges as well [5]. Therefore the limit m — > gives 
access to the jammed glassy states at infinite pressure [5]. 

The limit for m going to zero of Z™ nk , ZJ£ X , and Z™ te are either zero (for "incompatible" configurations 
of the cavity fields) or one (for "compatible" configurations of the cavity fields) regardless of the value of the 
cavity fields. As a consequence, in order to compute the complexity (which equals the replicated entropy 
S(m — > 0), since the internal entropy term, ms(m), disappears) we are only interested in the propagation of 
this information. 

To this aim, we introduce the "hard" components of the cavity fields iphard and Phard'- 

. , > f 1 if ip(x) > , , / 1 if <p(x) > , . 

Vw(z) = I otherwise and p hard {x) = | Q otherwise (78) 

These functions are defined as being equal to one for all values of x such that the cavity fields are non 
vanishing regardless of their value (i.e., corresponding to a non- vanishing probability of finding a sphere with 
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center in x), and zero otherwise. Since the reweighting factors in Eq. (13) do not depend on the actual value 
of the fields in the m — > limit, the propagation of the hard components decouples completely from the 
propagation of the cavity fields and can thus be treated indepenently. As a consequence, the population 
dynamics algorithm described above can be used on the populations encoding the probability distributoons 
of the hard fields. Once a stationary state has been reached, we can compute the complexity at m = 0, 
£j, from Eq. (12), computing the logarithm of the average value of the fraction of attempts yielding a non 
vanishing value of Zi in k, Zb ox , and Z site . Using Eq. (73), instead of computing (Z^ x ) and (Z™ te ), one can 
more easily compute (Z™) and {Z™), which are given respectively by the average value of the fraction of 
non-rejection attempts to generate the new tphard and iphard fields over the total number of attempts. Then 
we can determine the location of Dqcp defined as Y*j(Dccp) = 0. 

The results at m = obtained with this method will be reported in Sec. VIII, and compared with the 
analytical approximations. 

An important caveat is that in principle some fields could be proportional to exp(— 1/m) in the limit 
m — > 0. If this happens, then the procedure above fails since these fields give a finite contribution to the 
normalizations which is neither nor 1. Although we could not perform a careful systematic investigation 
of this effect, it seems that it might happen only for values of z and p where the transition at m = 1 is 
continuous. This point surely deserves further investigation. 

Note that in order to compute the correlation function in the close packing limit (see Sec. IX) we also 
need to know the actual values of the cavity fields. Since the propagation of the hard components decouples 
completely from the the one of the fields itself, one can use the population dynamics algorithm to find the 
solution of the 1RSB equations for the distributions of hard fields and of the cavity fields independently 
(knowing that the cavity fields can only be non zero where the hard components are equal to one), and use 
Eq. (80) to compute the pair correlation function. 

VIII. COMPARISON BETWEEN NUMERICAL RESULTS AND THE APPROXIMATIONS 

In this section we report the results obtained from the direct numerical calculation with discrctized space 
and we compare them with the delta and Gaussian approximations. 

A. Complexity 

In Fig. 5 we report the complexities Y> eq (the complexity at m = 1 equal to (l/N) time the logarithm of 
the typical number of glass states when configurations are samples uniformly) and Ej (the complexity at 
m = equal to (1/N) time the logarithm of the total number of jammed states) for several representative 
cases at d — 1 where the transition is discontinuous. Generically we observe that the delta approximation 
performs better at m = 0, while the Gaussian approximation is more reliable at m = 1. Both approximations 
give an upper bound to the true complexity and therefore give values for £>k and -Dgcp that are above the 
true ones. Moreover, both approximations miss the dynamical transition since by construction the fields are 
assumed to be localized. 

Some results for d = 2 are reported in Fig. 6. Here the scaling for q — > oo becomes very difficult 
because the numerical solution is computationally demanding and we cannot go beyond q = 20 for moderate 
connectivities. We could perform a systematic investigation only p = 2 and z = 20, which is unfortunately a 
case where the transition is continuous and the solution might be unstable towards further RSB in the glass 
phase. In this case, at m — 1 we correctly find a continuous transition at a value of D which is compatible 
with the result found from the stability analysis of section IV. At in = 0, we find good agreement with 
the result of the Gaussian and delta approximation. Note however that also at m — the results could be 
unstable towards further RSB. 
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FIG. 5: The complexity in some representative cases of discontinuous transition at d — 1, computed with the numerical 
solution of the population dynamics algorithm with varying resolution of the discretization process, is compared to 
the Gaussian and the Delta approximations. [Upper panels] E e9 (left) and Ej (right) for p — 2 and z = 110. In both 
cases we fixed the parameter a = 4, 7, 10 in Eq. (69) and changed q to vary the effective diameter D e g = (l + 2a)/(2g), 
which is reported in the horizontal axis. [Lower panes] E eg (left) and Ej (right) for p — 4 and 2 = 3. In the first 
case, we varied q at fixed a, while in the second we did the inverse. 



B. Phase diagram 



In Fig. 7 we compare the transition lines obtained by the Gaussian and delta approximations with the 
numerical results, where available. We computed _Dk and -Dgcp by performing an extrapolation to q — > oo 
(which is simple since the corrections are found to be proportional to 1 /q) in some representative cases where 
the transition is continuous or discontinuous; the results are reported in Fig. 7. We observe that indeed the 
Gaussian and delta approximation give consistent results, which are also consistent with the exact numerical 
solution and provide upper bounds to the latter. 

Whenever the RS instability Z?r.s < £*k, the transition is continuous. This happens generically for small 
z. On increasing z, the lines Drs and Z?k cross and the transition becomes discountinuos. The value z* 
where this crossover happens depends weakly on the space dimension, but it depends strongly on p. Indeed 
we have z* ~ 100 for p = 2, while z* ~ 20 for p = 3 and (as we can infer from Fig. 4) the transition is always 
discontinuous for p > 3. 



25 




FIG. 6: The complexity at d = 2, p = 2 and « = 20, computed with the numerical solution of the population 
dynamics algorithm with varying resolution of the discretization process, is compared to the Gaussian and the Delta 
approximations. Here we can only use moderate values of q, and because of the geometry of the discretization the 
effective diameter of the sphere, given by Eq. (70) and reported on the horizontal axis, cannot be varied smoothly. 
For instance, at q = 11 we could not find a point at positive complexity. 



IX. CORRELATION FUNCTION 



A. Definition 



As explained in section III, in the glass phase the cavity equations have multiple solutions, each describing 
a different glass state. Within each state a we can define a correlation function g a {x, y) as follows. For each 
box we have: 



(79) 

1 J dx\ ■ ■ ■ dx a p Vff K)- (ag) • • • , x a p ) £ S(x - xf)S(y - xj) 
PiP- 1 ) Jdxt---dx-^l(xt)---^(x-)x(xt,--- ,z») 

since the fields ip^i(xf) describe the distribution of the variables adjacents to box a in absence of the box 
itself. We now average this quantity over the boxes and over the states a with the weight Z™ . We get 

Nz/p 

JVZ =1 ^a^a a 

(80) 



J dV[ipi}---dT[ipp}Z box [ipi---'ipp} m 1 ipi(x)ip2(y) J ^f[^j{xj)dxjjx(x,y,X3,---,Xp) 
Note that in the RS case the above expression reduces to g~(x, y). 



K j=3 




FIG. 7: Phase diagrams for p = 2, 3 and d = 1, 2. We compare the results of the Gaussian and delta approximations 
with the numerical results obtained directly from a discretization of the cavity equations. In the lower right panel, 
the horizontal line indicates the value D = 1/4 above which the calculation of Z$ is not valid, see Eq. (29). 



We expect that at m — (close packing), g(x, y) develops a peak in \x—y\ = D describing contacts [27, 28]. 
The number of contacts is 

C = (p-i)/ 9(0,y)dy. (81) 

J peak 

The delta peak is also accompanied, in three dimensional sphere packings, by a square root divergence, 
g(r) ~ (r — D)~ 5 [27, 28], which we want to investigate here. 
Note that in the delta approximation we just get 

9(x, y) = -^o J dX 3 ■ ■ ■ dx pX(x, V,X 3 ,--- , X p ) = g°(x, y) (82) 

therefore all the structure of the correlation in the packings is lost in this approximation. 

One can show, following [5], that in the Gaussian approximation, as A ~ m for m — > 0, one gets a 
delta peak at r = D in the jamming limit, with all particles being non-rattlers and £ = 2d. Therefore this 
approximation is able to capture some of the peculiar structure of the correlation. On the other hand, the 
square root singularity is missed by the Gaussian approximation [5]. 

Unfortunately, it is very difficult to study the contact peak in the numerical solution of the cavity equation, 
because the discretization makes it hard to define a proper notion of contacts and separate the delta peak 
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FIG. 8: Pair correlation function g(r) at d = 1, m = (jamming) and D ~ -Dgcp (in practice, the closest value 
to -Dgcp compatible with the discretization). (Left) p — 2, 2 = 6; note that in this case the system undergoes a 
continuous transition and these results might be unstable towards further RSB. (Right) p = 4, 2 = 3: here the 
transition is discontinuous. Note that for p = 4 we observe an additional singularity at r — 2D [27]. 



contribution from the background. Therefore in the following we focus on the square root singularity which 
is also a non-trivial and somehow unexpected feature of pair correlations at jamming [27, 28]. 

Numerical results are presented in Fig. 8 for the g(r) in one dimension, and two representative values of z 
and p where the transition is continuous or discontinuous. In both cases, the divergence is compatible with 
a square root singularity (r — _D)~ 0,5 in a range of r — D, but at smaller r — D the g(r) seems to diverge as 
(r — D)~ J with an exponent 7 > 0.5. However, in this region the square root divergence is probably mixed 
with the contact delta peak, because of the discretization. A detailed analysis of this mixing was not possible 
because the values of q we could reach were still too small. Since this investigation is computationally very 
demanding, we could not perform a systematic study of the value of the exponent as a function of p and 
z, nor investigate the more interesting case d — 2, which is very hard because our discretization does not 
preserve the spherical symmetry around the central particle. We leave a more systematic numerical analysis 
for future work. 



B. Argument for the square-root singularity 
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We now present an analytical argument to relate the shape of the cavity fields to the square root singularity. 
We focus on m = 0, and we study the small r — D behavior of g(r) as follows. We define the quantity 

*(*) = / (f[M^)dx)j X{x ^ x j* p) 6( xi -x 2 -z) . (83) 

Note that z = x — y g [— 1, l] d but using periodicity one can restrict to z g [— l/2,l/2] d with periodic 
boundary conditions. The probability distribution of tp induces a distribution V[&] on vp. Then we have 



g(z) = J dxdyg(x, y)S(x - y - z) = e "*~ J dV[*} r^^\j J dz^(z) X (z) 



(84) 



where the term e~ Sbo:L ensures the normalization J dzg(z) = 1. 

In the following we restrict for simplicity to d = 1. Note that by translational invariance the field \I> 
is centered around a random uniformly distributed position z , while its shape is encoded by a non-trivial 
distribution. Now assume that with a certain finite probability with respect to the shape distribution, one 
has that 

• ty(z) vanishes at some finite distance from the center given by z± — zq ± 5zq. The quantities z± are 
then also random and uniformly distributed in [—1/2, 1/2]; 

• the shape of &(z) around the point where it vanishes is of the form 



- e !"± |0 ; (85) 

• and \z+ — z- \ < 2D, and z + > D (the additional symmetric contribution coming from z_ gives a factor 
2 and will be neglected as all proportionality constants). 

Then the function x(z)^(z) vanishes everywhere except in [D, z+] where it is given by exp [ — A/ (z + — z) Q ] . 
The average over V[^} 7 for what concerns this contribution, is translated onto an average over z + and 
Eq. (84) becomes : 

r e~~ A / (z +~ z ^ a 0(D < z < z+) f c e -A/(z+-z) a 

» W ~/^ J^e-.UT-r ^ Zt-J. ^ J^m^F ■ < 86 > 

where C is a suitable cutoff that comes from the fact that if z + is too much larger than D the approximation 
Eq. (85) will break down. We will show that this cutoff does not matter as the main contribution for z — >• D 
comes from z + close to D. 

To simplify notations, we introduce A = (z — D)/D and e = (z + — D)/D. Also we define a = A/D a and 
c = (C — D)/D. With these notations we get 

r c e -a/(e-\) a 
5(A) ^ A d£ ^dXe-^r • (87) 

The integral in the denominator is dominated by the small A behavior, that gives 

f d\e- a ' {£ - x ^ - /" dXe-^( 1+a ^) = e -^ 1 —- , (88) 
Jo Jo aa 

and 

3 (A)oc / dee- {a+1 ^e a i^-T^w) . (89) 
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We want now to evaluate the integral by a saddle point for A — > 0. We assume (and will check self-consistently) 
that the saddle point value e* » A. Then we can expand for A/e <C 1 and 

g(X) oc J" de e- (a+ V ^s^-^- (a+1 ' > . ( 90 ) 

The maximum of the above expression is found at e* — (aaA) 1 ^"^ 1 -* 3> A for small A as initially assumed. 
Substituting this in the expression above one obtains 5(A) oc 1/A. To get the correct result we need to 
compute also the quadratic corrections around the saddle point. Including these, we finally obtain 

g(X) oc A"Tfe a (r - D)~^ , (91) 

i.e. a power-law divergence for z — > D with exponent € [0, 1], which is consistent with the observed exponents 
in Fig. 8. Note that a square root singularity is obtained for a = 1, namely a simple exponential singularity 
of the cavity fields. We checked on our numerical results that indeed the form of the fields is compatible 
with the Ansatz (85). 

Note that this same argument can be carried out at finite m, but in this case we get that g(X) is independent 
of A for small A. A more complete analysis should show that at finite m, g(X) is a power law for A » 0{m v ) 
with some exponent v, and it crosses over to a finite value for A <C 0(m' / ). 

X. DISCUSSION ON FINITE DIMENSIONAL HARD SPHERES 

One way to recover the normal hard sphere model from our model is to set p = 2 and z — N — 1. However, 
this limit cannot be investigated within the cavity formalism which is based on taking first the limit N — > 00 
at finite z. Here the limits N — > 00 and z — > 00 do not commute, and if we first send N — > 00 and then 
z — > 00 we do not recover the hard sphere models (a similar behavior is found for the Bethe lattice spin 
glass [23]). 

Therefore we want here to find a suitable limit that we can take after N — > 00 to recover the hard sphere 
model. As we discussed in the introduction, one possibility if to set formally z = 1 and identify p with the 
number of particles, therefore taking p> 1. Of course, for z < 2 and finite p the model does not have any 
phase transition (it becomes a one-dimensional model for z = 2). Therefore, we have to send p — > 00 before 
z becomes smaller than 2. 

As a first check, we note that in this limit the RS entropy 

S RS = Z -\ogZl^SuM , (92) 

where Sn q (ip) is the entropy of d-dimensional hard spheres in the thermodynamic limit at fixed packing 
fraction cp. Actually, there is a problem with the latter identification, since Z® does not contain a factor 
pi which should take into account indistinguishability of the particles. This is indeed to be expected, since 
we took a formal limit z — > 1, but at any finite z > 1 the particles are connected to several boxes which 
makes them distinguishable. We therefore recover the finite dimensional result for a system of distinguishable 
particles. 

Next, we can look at the stability of the RS solution according to Eq. (20). To compare with standard 
hard spheres it is crucial to observe that here the box side is one while D becomes very small for p — > 00, 
in such a way that the packing fraction ip = pVd(D/2) = pVd(l/2)D d is finite. For p — > 00 first and z — > 1 
after, we have <7p(x) — > gu q (x), however x is expressed in units of the box length. If we introduce as usual 
the distance r measured in units of the sphere diameter, r = x/D, we have (for k ^ 0) 

g° p (k) = J dxe lkx g° p (x) = D d J dre lkDr g liq (r) = D d S{kD) , (93) 

where S(kD) is the structure factor, and the stability condition becomes 



^-1)^-1)^15(^)1 = ^(P -!)(*-!) * lo A S(kD)\ < 1 



(94) 
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which is always verified for p — > oo since ip and S(kD) are both of order 1. This is indeed consistent with our 
investigations of the model at finite p that showed that the transition is always discontinuous at p > 4. We 
conclude that one cannot observe a continuous transition in the normal hard spheres model. This conclusion 
is consistent with the ones of Biroli and Bouchaud [29] who showed that indeed replicated liquid theory in 
finite dimensions does not allow for a continuous RSB transition. 

We also note that starting from Eqs. (43), (42) and taking first p — > oo (with <p = pVd(D)/2 d and 
(p— l)<7p(r) = pgu q {r)) and then z^lwe recover Eq. (74) of [5], which is the starting point of the Gaussian 
small cage replica treatment in finite dimensions, provided we identify again lim p _>. 00 MogZ^ = Su q (ip), 
neglecting the problem with the missing pi. Apart from this caveat, this is a nice alternative derivation of 
the approximation of [5], which is not based on the replica method. 

Finally, one could try to take the same formal limit in Eq. (52) to obtain an alternative approximate 
expression for S(m) in finite dimensions. Using the relation Zp/Zp_ 1 = (v), where v is the void space of 
p — 1 particles, we obtain for z — > 1 (after p — > oo): 

( v m \ 1 

S(m) = IogV-r + -Iog^. (95) 
(v) p 

Note however that the void space v oc p, therefore we must rearrange terms as 

S(m) = l0g +ml ° gP + \ Xo ^ Z pIp P ) ■ (96) 

The term mlogp can be dropped since it gives an additive constant to the internal entropy, and the resulting 
expression has a well defined p — >• oo limit, assuming here that linip^oo ilog(Zp/p!) = Su q (<p) (which is 

however inconsistent with the previous discussion, for reasons that we do not understand at present). This 
expression can in principle be directly computed, even if it is very hard to sample the distribution P(v) of 
void space because at high density v = for most configurations [30]. 



XI. CONCLUSIONS 



In this paper, we have studied a mean field hard sphere model introduced in [12]. The model is similar 
to a standard hard sphere model, however each sphere interacts only with a finite and preassigned number 
of neighbors. The network of interactions is given by a random graph, such that the model belongs to the 
mean field class and is therefore, in principle, exactly solvable via the cavity method. We therefore derived 
the cavity equations for the model and we presented both analytical approximations to their solution and 
an "exact" numerical solution based on a discretization of the space. 

We have shown that the analytical approximations give quite reliable results for the phase diagram and 
the complexity. In particular, for large enough z and/or p, the transition belongs to the Random First Order 
class. Therefore, as suggested in [12], the model displays an ideal glass (Kauzmann) transition to a glass 
phase. Following the glass phase upon increasing pressure, one gets to a point where the pressure diverges, 
similarly to standard hard spheres close to the so-called J-point. Given that the model has an exponential 
number of metastable states, one obtains a set of J-points spanning a finite range in density. Overall, the 
phenomenology of the model in this regime is very close to the one expected for finite dimensional hard 
spheres based on mean field approximations, see [5] and Fig. 1. We found, in particular, that the Gaussian 
approximation is very good for the Kauzmann transition but tends to overestimate the close packing. This 
is consistent with what happens for three-dimensional hard spheres where the Gaussian approximation gives 
(fK ~ 0.62, which is consistent with numerical estimates, and (fccp ~ 0.68, while numerical simulations 
suggest a somewhat smaller value [5]. On the contrary, the delta approximation is very good for close 
packing but tends to overestimate the Kauzmann point. We proposed a formula for the complexity that is 
based on the delta approximation and can be computed numerically for three-dimensional hard spheres. It 
would be very interesting to do this computation and compare the result with the Gaussian approximation 
in that case. 
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We also found a somehow unexpected result, that the transition is continuous at small z and p. In 
particular, for the values of p = 2 and z = 100 that have been used in [12], the transition should be very 
weakly first order. The physics in presence of a second order transition could be very different. For instance, 
in the case of the Sherrington-Kirkpatrick model, the intensive ground state energy can be found easily: this 
would correspond to a unique J-point density. However, the details of this depend on the model, and in 
particular on the shape of the complexity function, so we cannot give any conclusive statement. It would be 
interesting to investigate better this point by repeating the numerical simulations of [12] both in a region 
where the transition should be strongly second order (e.g. at p = 2 and small z) and in a region where it 
should be strongly "random first order" (e.g. for p = 4 and small z). 

Finally, we partially investigated the structure of the configurations at jamming. We computed the 
correlation function of the model and showed that it displays a power-law singularity close to contact, 
at least for d = 1 . We also gave an analytical argument to explain the mathematical origin of the singularity. 
Extending this study to higher dimension could give insight in the physics that is responsible for this 
divergence and hopefully connect it to isostaticity and the presence of soft modes in the spectrum, as 
suggested in [14, 15]. Additional numerical simulations could be extremely useful also in this respect. 
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